function [log_likelihood, nan_pct, sim_factors]=ml_separate_factors(theta)
%%% calculates the log-likelihood allowing multiple factors. 

 global f_ijs X_ind_s measures ind_id factor_n  ind_n draw_n c_char_n b_char_n m_char_n cont_char_n thresh X_n theta_xy lottery_n lottery_choice_ind time_choice_n f_n
    
    %%% determines factor being estimated
    i=f_n;

    % each measure is associated with a unique mean and loading coefficient
    % for the factor   
    f_measures=measures(:,sum(c_char_n(1:i-1))+1:sum(c_char_n(1:i-1))+c_char_n(i));
    columns=sum(c_char_n(i));
    % # of observations including rows with simulated factor draws for each
    % person
    rows=size(f_measures,1);
    
    sex=X_ind_s(:,1);
    
        
    % gets coeffs on the demographic vars for each factor 
    alphas=theta(1:X_n);
    theta(1:X_n)=[];

    % all multi-valued variables have the same amount of categories
    % for a given factor
    % number of nuisance parameters (thresh-mean) to be estimated for
    % this particular factor
    size_const=b_char_n(i)+thresh(i)*m_char_n(i)+cont_char_n(i);
    %%% intercept for 1st measure of each factor is normalized to 0
    constants=[0; theta(1:size_const-1)];
    %%% factor-specific intercept
    f_mean=theta(size_const);
    % one loading per measure
    loading=[1 theta(size_const+1:size_const+c_char_n(i)-1)'];
    %%% factor orthogonal component sd
    sd=exp(theta(size_const+c_char_n(i)));
    %%% error sd for the continous measures
    sd_cont=exp(theta(size_const+c_char_n(i)+1:size_const+c_char_n(i)+cont_char_n(i)));
    theta(1:size_const+c_char_n(i)+cont_char_n(i))=[];        
    %%% to be filled in gradually with the binary, multi, and
    %%% continuous measure likelihood contributions
    sim_l_contrib=zeros(rows,c_char_n(i));

    %%% "f_ij" are simulated factor draws (common to all measures)
    %%% associated with a given factor
    sim_factor=f_mean+X_ind_s*alphas+f_ijs(:,i)*sd;


    if b_char_n(i) > 0

        % gets the binary measures for the factor
        X_c_b=f_measures(:,1:b_char_n(i));
        %nuisance parameters for the binary measures associated with the
        %factor
        t_k_mean_b=constants(1:b_char_n(i))';
        % remove the already assigned parameters from the constant list
        constants(1:b_char_n(i))=[];
        % loadings for binary measure
        loading_b=loading(1:b_char_n(i));
        % remove the already assigned parameters from the loading list
        loading(1:b_char_n(i))=[];



        %%%%%%%%%%%% Simulated likelihood contributions of the BINARY
        %%%%%%%%%%%% measures 

        %%% simulated measure value for each individual and simulated factor
        %%% draw - excludes the constant term for multi-valued categorical
        %%% estimation
        sim_u=sim_factor*loading_b;
 
        temp=bsxfun(@plus,-sim_u,t_k_mean_b);       
        sim_l_contrib(:,1:b_char_n(i))=((1-normcdf(temp)).^(X_c_b)).*(normcdf(temp).^(1-X_c_b));

    end






    %%%%%%%%%%%% Simulated likelihood contributions of the MULTI-VALUED
    %%%%%%%%%%%% measures (calculates sim_u and sim_l_contrib for those now)

    if m_char_n(i) > 0
        %%% gets the multi-valued measures for the factor
        X_c_m=f_measures(:,b_char_n(i)+1:b_char_n(i)+m_char_n(i));

        t_k_mean_raw=zeros(thresh(i),m_char_n(i));
        for j=1:thresh(i)
            t_k_mean_raw(j,:)=constants(1+(j-1)*m_char_n(i):j*m_char_n(i));
        end

        t_k_mean_m=zeros(thresh(i),m_char_n(i));
        t_k_mean_m(1,:)=t_k_mean_raw(1,:);
        % ensures monotonicity in threshold-mean coeffs for each measure at all
        % times (not just at start) 
        for j=2:thresh(i)
            t_k_mean_m(j,:)=t_k_mean_m(j-1,:)+exp(t_k_mean_raw(j,:));
        end  

        % remove the already assigned parameters from the constant list
        constants(1:thresh(i)*m_char_n(i))=[];

        % loadings for multi measure
        loading_m=loading(1:m_char_n(i));
        % remove the already assigned parameters from the loading list
        loading(1:m_char_n(i))=[];

        %%% simulated measure value for each individual and simulated factor
        %%% draw - excludes the constant term for multi-valued categorical
        %%% estimation
        %%% "f_ij" are simulated factor draws (common to all measures)
        %%% associated with a given factor
        sim_u=sim_factor*loading_m;

        % probability contribution of multi-valued measures: 
        % lower end
        sim_l_contrib_m=zeros(size(X_c_m,1),size(X_c_m,2));
        temp=bsxfun(@plus,-sim_u,t_k_mean_m(1,:));
        sim_l_contrib_m(X_c_m==0)=normcdf(temp(X_c_m==0));
        %thresholds in between
       for j=1:thresh(i)-1
            temp=bsxfun(@plus,-sim_u,t_k_mean_m(j,:));  
            temp1=bsxfun(@plus,-sim_u,t_k_mean_m(j+1,:));         
            sim_l_contrib_m(X_c_m==j)=normcdf(temp1(X_c_m==j))-normcdf(temp(X_c_m==j));
       end
       % upper end
       temp=bsxfun(@plus,-sim_u,t_k_mean_m(thresh(i),:));  
       sim_l_contrib_m(X_c_m==thresh(i))=1-normcdf(temp(X_c_m==thresh(i)));

       sim_l_contrib_m(isnan(X_c_m))=1;

       %%% Adding the multi-valued contributions to the overall factor likelihood 
       sim_l_contrib(:,b_char_n(i)+1:b_char_n(i)+m_char_n(i))=sim_l_contrib_m;

    end


    %%%%%%%%%%%% Simulated likelihood contributions of the CONTINUOUS
    %%%%%%%%%%%% measures (calculates sim_u and sim_l_contrib for those now)
    %nuisance parameters for the multi-valued measures associated with the
    %factor
    if cont_char_n(i) > 0

        %%% Gets the continuous measures for the factor 
        X_c_c=f_measures(:,b_char_n(i)+m_char_n(i)+1:b_char_n(i)+m_char_n(i)+cont_char_n(i));
        %nuisance parameters for the continuous measures associated with the
        %factor
        mean_c=constants(1:cont_char_n(i))';
        % remove the already assigned parameters from the constant list
        constants(1:cont_char_n(i))=[];

        % loadings for continuos measure
        loading_c=loading(1:cont_char_n(i));
        % remove the already assigned parameters from the loading list
        loading(1:cont_char_n(i))=[];        

        % transforms the mean vector to a matrix so that it can be added to he
        % simulated loading contribution (the mean is common to all persons for
        % a given measure)
        mean_f=repmat(mean_c, rows,1) ;


        % simulated factor draws (common to all measures)
        sim_u=mean_f+sim_factor*loading_c;

    
        %%% "mu" just says that the error terms of each measure are assumed
        %%% mean zero (and sd sigma but that comes in later)
        mu=zeros(size(sim_u));
        %%% allows for a measure-specific error term 
        sigma=repmat(sd_cont',rows,1);
        sim_l_contrib(:,b_char_n(i)+m_char_n(i)+1:c_char_n(i)) = pdf('Normal',X_c_c-sim_u,mu,sigma);       
    end    
        

  
    
   %%%%%%%%%%%%%%%%%%%%%% Combining the contributions from the various
   %%%%%%%%%%%%%%%%%%%%%% types of measures into a likelihood
   %%%%%%%%%%%%%%%%%%%%%% contribution for each individual and factor
   %%%%%%%%%%%%%%%%%%%%%% draw
    
   sim_l_contrib_total=[sim_l_contrib];
   %%% for numerical optimization purposes
   sim_l_contrib_total_plus=sim_l_contrib_total+1*10^(-20);
   
    % Individual likelihood contribution conditional on a factor draw
    sim_ind_draw_contrib=prod(sim_l_contrib_total_plus,2);

    % individual average of product of measures per draw and lottery
    % choices
    av_l_contrib=cell2mat(accumarray(ind_id,sim_ind_draw_contrib,[],@(x){sum(x)}))/draw_n;

    % In order to avoid zeros for the log - enough that one of the
    % likelihoods in a row is zero and whole row is out otherwise and
    % causing trouble for the log (can easily happen as the minimization
    % routine goes through various parameters, notably with random starting
    % values)
    av_l=av_l_contrib+1*10^(-200);    
    
    log_likelihood=-nanmean(log(av_l));
    
    %%% will capture the % of nan and inf in the final average likelihoods
    nan_pct=1-mean(isfinite(av_l_contrib));
    if nan_pct>.05
        log_likelihood=NaN;
    end
    
end